Nonequilibrium Green's functions approach to strongly correlated 

few-electron quantum dots 

K. Balzer and M. Bonitz 

Institut fiir Theoretische Physik und Astrophysik, 
Christian- Albrechts-Universitat Kiel, Leibnizstrasse 15, 24098 Kiel, German^ 

R. van Leeuwen 

Department of Physics, University of Jyvaskyla, Survontie 9, 4OOI4 Jyvdskyld, Finland 

N.E. Dahlen and A. Stan 
Zernike Institute of Advanced Materials, University of Groningen, 
Nijenborgh 4, 9747 AG Groningen, The Netherlands 
(Dated: October 14, 2008) 

The effect of electron-electron scattering on the equilibrium properties of few-electron quantum 
dots is investigated by means of nonequihbrium Green's functions theory. The ground and equi- 
librium state is self-consistently computed from the Matsubara (imaginary time) Green's function 
for the spatially inhomogeneous quantum dot system whose constituent charge carriers are treated 
as spin-polarized. To include correlations, the Dyson equation is solved, starting from a Hartree- 
Fock reference state, within a conserving (second order) self-energy approximation where direct and 
exchange contributions to the electron-electron interaction are included on the same footing. 

We present results for the zero and finite temperature charge carrier density, the orbital-resolved 
distribution functions and the self-consistent total energies and spectral functions for isotropic, two- 
dimensional parabolic confinement as well as for the limit of large anisotropy — quasi-one-dimensional 
entrapment. For the considered quantum dots with A'^ = 2, 3 and 6 electrons, the analysis comprises 
the crossover from Fermi gas/liquid (at large carrier density) to Wigner molecule or crystal behavior 
(in the low-density limit). 



PACS numbers: 73.21.La, 05.30.-d, 71.27.-|-a 
I. INTRODUCTION 

In the recent decades finite quantum systems have be- 
come an intensively studied subject-matter. Particular 
interest is due to electrons in quantum dots'^ (QD) or 
wells, forming artificial atoms^ with molecule-like behav- 
ior and novel spectral and dynamical properties. In con- 
trast to real atoms, such new properties arise from dimen- 
sionality reduction and natural scale-differences as QDs 
embedded within semiconductor heterostructures gener- 
ate charge carrier motion on typically nanometer length 
scales. The collective, optical and transport properti es of 
QDs are examined by experimental^ and theoreticaP'^ 
research activities in dependence on various dot p aram- 
eters and geometries. For an overview see e.g. Refs.l^^^. 
Many ground state calculations are available in the lit- 
erature which are based on d ifferent methods — exact 
num erical d iagonalizatioij ^ l ^° l ^^ l, self-consistent Hartree- 
PQpj^li2li3ii4j^ configuration interactiorfl^ and quantum 
Monte CarlcJi^. Extensions to finite temperatures and 
to QD propertie s in (tr ansversal) magnetic fields are to 
be found in Refs! ^ ' ^^ ' " l. 

Typical charge densities in QD devices vary over a 
large rang^^^l — from macroscopic charge carrier ensem- 
bles to mesoscopic, few- and even single-electron quan- 
tum dots. However, injecting only a small integral num- 
ber of electrons into the dot reveals system properties 
that sensitively depend on the charge carrier number 



and are thus externally controllable, e.g. by gate voltage 
or tip-electrode field variation or local mechanical strain 
(band mismatch). On the other hand, the quantum dot 
state is governed by the interplay of quantum and spin ef- 
fects, the Coulomb repulsion between the carriers and the 
strength of the dot confinement. This generally leads to 
strong electron-electron correlation i.e. collision or scat- 
tering effects the influence of which on the many-particle 
state is very important in the behavior at zero and finite 
temperature. 

The two-dimensional (2D) A^-electron quantum dot 
Hamiltonian to be considered is 

where the effective electron mass is denoted by m*, 
the frequency wq adjusts the (isotropic) parabolic con- 
finement strength, e is the elementary charge and e is 
the background dielectric constant. The vectors are 
the single charge carrier coordinates with respect to the 
quantum dot center and rij = jr^ — rj|. The density in 
the QD is controllable by the confining potential which 
directly affects the relative electron-electron interaction 
strength and, tuned towards low carrier densities, contin- 
uously le ads to formation of so-called electron (Wigner) 
molecule^i^^ or crystal-like behavioi^^. Melting pro- 
cesses owing to an increased temperature cause weaken- 
ing and finally preventing of such solid-like structure for- 
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mation. Also, with the restriction to Eq. (jlj, the present 
analysis neglects nonideality effects such as defects and 
well-width fluctuations, see e.g. 

The objective of the present work is to analyze the 
electron correlations in the QD system ([!]) on the path- 
way from Fermi gas/liquid towards strongly correlated 
Wigner molecule behavior, i.e. during the (density) 
delocalization-localization transition. To this end, in 
Sec. [n] we present the finite temperature formalism of 
nonequilibrium Green's functions (NEGF) theory which 
is applied by a threefold motivation: (i) the NEGFs al- 
low for a consistent and conserving treatment of Coulomb 
correlations, (ii) previous NEGF approaches to inhomo- 
geneous QDs incorporate to our knowledge no strong 
carrier-carrier coupling, and (iii), in contrast to other 
methods, strong importance lies in the possibility for the 
direct extension of the approach to nonequilibrium situ- 
ations with e.g. time-dependent gate-voltage variations, 
quantum transport phenomena or optical switchinj^l. 

Using the NEGF technique, the properties of the inves- 
tigated spin-polarized A^-electron quantum dot in ther- 
modynamic equilibrium follow from the self-consistently 
obtained imaginary time (Matsubara) Green's func- 
tion. Such an approach has also more recently shown 
to give accurate results for real atoms and molecules, 
see Refs.l^fllSllSlI ^j^g extension of the nonequilibrium 
Green's function ansatz from traditi onal app lications on 
quasi-homogeneous quantum systemg^^ISiESlto spatial in- 
homogeneity is thereby the major goal of the present 
analysis. 

The theoretical part in Sec.|ll]is followed by a detailed 
description of the iteration technique used to numerically 
solve the Dyson equation (in Hartree-Fock and second 
Born approximation) to self-consistency. The results are 
discussed in Sec. |IV[ The starting point is the limit of 
large anisotropy where in Eq. ([T]) we consider the limit 



II. THEORY 



ujqV^ ^iflX^ + ^yflVi witii ijjyfl > cja;,o, see Sec. IV A 
In the case of = 3 (quantum dot lithium) and d eiec- 
trons, the charge carrier density, the orbital-resolved dis- 
tribution functions and total energies are computed for 
different values of interaction strength and temperature. 
Also, we demonstrate that at finite temperatures, the 
second Born (correlation) corrections to the mean-field 
treatment yield significant density changes in an inter- 
mediate regime whereas in the high- and low tempera- 
ture limit the electron density is only less affected by 
correlations. In Sec. |IVB[ we extend the calculations 
to isotropic 2D confinement and analogously report on 
ground state results for N — 2 electrons (quantum dot 
helium) which are compared with exact and Hartree- 
Fock results of Ref.^^. Moreover, the computation of the 
charge carrier spectral functioiP^a(a;) allows in Sec. IV C 
for a collision induced renormalization of the Hartree- 
Fock energy spectrum. This is of high relevance for the 
optical properties of the few-electron QD. Sec. |V] gives a 
final discussion. 



For characterization and quantum mechanical treat- 
ment of the A-electron dot system ([T|) it is convenient to 
introduce the coupling (or Wigner) parameter A which re- 
lates the characteristic Coulomb energy Ec = e^/(47re/(5) 
to the confinement energy Eq = Hloq : 



Er 



1* 



Eq iTTelQhujo as 



(2) 



with Iq 



yjh/ {m*ujQ) being the characteristic single- 
electron extension in the QD and as the effective elec- 
tron Bohr radius. Using the replacement rules {r^ 
Yi/lQ, E E/Eq}, Hamiltonian (jlj transforms into the 
dimensionless form 



1 ^ 

i=l 



N 

E 



(3) 



For coupling parameters A ^ 1, the quantum dot elec- 
trons will be found in a Fermi gas- or liquid-like state, 
whereas in the limit A — s- oo, it is Zq S> as, and quantum 
effects vanish in favor of classical interaction dominated 
charge carrier behavioil^. In the case of moderate cou- 
pling (A > 1) quantum dots with spatially well localized 
carrier density can be formed. Further, in addition to N 
and A, the system is characterized by the QD tempera- 
ture (3~^ = ksT which will be measured in units of the 
confinement energy Eq. 



A. Second quantization representation 

Introducing carrier annihilation (creation) operators 
'0(t)(]-) with action at space point r, the second-quantized 
form of Hamiltonian H\, Eq. is 



H)^ = /d^rV;^(r)/i°(r)V;(r) 



+ i lldV d^f (r) (r) ^=|p 



(4) 

-0(f) ■0(r) , 



where h^{Y) — (— + r^)/2 denotes the single- 
particle energy and the second term in ^ describes 
the electron-electron interactions. The field operators 
'^^^•'(r) satisfy the fermionic anti-commutation relations 
[^/'(r),V'^(f)]^ = 5{v-v) and [V>(t)(r)^ ^(t)(f )]^ ^ q, 

where [i, 13]+ ^ AI3 + BA. 

Ensemble averaging in Eq. ^ directly gives rise to 
the one-particle nonequilibrium Green's function which 
is defined as 



G{l,2) = ~-{Tcm)ij\2)] 



(5) 



and is a generalization of the one-particle density matrix 
[which is recovered from G in the limit of equal time ar- 
guments ti = t2, see e.g. Ref.*^. The used nomenclature 
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is 1 = (ri,ti) and the expectation value (ensemble aver- 
age) reads {A) = Tr pA. The two times ti and t2 entering 
G{1,2) arise in the Heisenberg picture of the field op- 
erators and vary along the complex Schwinger/Keldysh 
contour C = {t e C\9\tt G [0, oo] , i e [-/3, Oil 
where Tq denotes time-ordering on C, see e.g. Refs.^^iJ 
ancpSl. Note, that in the remainder of this work we use 
h= 1. The advantage of using the NEGF is that it allows 
for equal access to equilibrium and nonequilibrium aver- 
ages at finite temperatures and that quantum many-body 
approximations can be syst emat ical ly i ncluded by dia- 
gram expansiond^, see Sees. II B and III Moreover, most 
dynamic (spectral) and thermodynamic informatiorP^ is 
contained in the NEGF, cf. Sec. llVl 



B. Nonequilibrium Green's functions formalism 

The two-time nonequilibrium Green's function G (l,2) 
obeys the Keldysh/Kadanoff-Baym equation (KBEpSMl 

[za,,-/j"(ri)]G(l,2) 
= Sc{l~ 2) d3 W{1 - 3) Gi2(13; 23+) , (6) 

and its adjoint equation [with interchanged time argu- 
ments ti t2]- On the right hand side of Eq. ^ 
the (collision) integral runs over the full configuration 
space and the time domain spanned by the contour 
C. Further, 1^(1 - 2) = XScih - t2) / {ri - r2)^ is the 
instantaneous (time-local) electron-electron interaction 
and (5c (1 ~ 2) = dc{ti — ^2) <5(ri — r2) with the time delta 
function being defined on the contour. Gi2(12; 1'2') de- 
notes the two-particle NEGF 

Gi2(12;l'2') = (-*)2(Tc[VXl)V'(2)V'^(2')V'^(l')]) , (7) 

where the short notation 3+ in Eq. ^ indicates that 
the limit t + is taken from above on the con- 

tour. In the integro-differential form ([6|, the KBE is not 
closed but constitutes the first equation of the Martin- 
Schwinger (MS) hierarchjP3. In order to decouple the 
hierarchy approximate expressions for the two-particle 
Green's function are introduced. E.g. in a first order 
(spatially non-local) Hartree-Fock approach one substi- 
tutes Gi2(12;l'2') -> G(1,1')G(2,2') - G(l, 2')G(2, 1') 
which is known to preserve total energy and momentum 
but completely neglects correlations — the former term 
leads over to the Hartree potential, the latter accounts 
for exchange. More generally, such conserving approxi- 
mations can be formulated in terms of a self-energy func- 
tional 5][G](1, 2) which is defined by 

-i J d3iy(l - 2)Gi2(13;23+) 
= ^d3S[G](l,3)G(3,2) . (8) 



Other, advanced conserving approximations, such as the 
second Born approximation (see Sec. Ill), can be sys- 



tematically derived from a generating functional <&[G] ac- 
cording to S(l,2) = 5$[G]/(5G(2,1), see e.g. Ref.^s. in 
addition to a specific MS hierarchy decoupling, the KBE 
^ must be supplied with initial or boundary conditions. 
In this p aper, we will use the Kubo-Martin-Schwinger 
conditioniPSl^ G(riti, 2)|t^=o = -G(ri - 2) and 
G(l,r2i2)|t.=o = -G(l,r2 0-1/3). 

In the case of thermodynamic equilibrium, where with- 
out loss of generality the electron system (jlj is time- 
independent for 9^eti_2 < 0, G(l,2) has no real-time de- 
pendence but extends on the imaginary contour branch 
[— z/3, 0] only. We define the corresponding Matsubara 
(imaginary time) Green's function G*^ with respect to 
the transformations ti — t2 — * ir (r e [— /3, /?]) and 



G 



~iG 



M 



G^>l,r2;T) = -^G(l,2) 



(9) 



which only depends on the time difference ti — t2, ^1,2 G 
[— i/3, 0] and is anti-periodic in the inverse temperature /3, 
compare with definition Using expressions ([8| and 
([9]) in the KBE (6]) leads to the general form of the Dyson 



equation^ for t 



re spin-polarized QD system M 



(10) 



[-a,-/i"'(ri)]G>i,r2;T) 
6{T)+Jd'fJ^ drEf(ri,f;T-f)G^^(f,r2;f) 



with the anti-periodic Matusbara self-energy 
I]*^(ri,r2; r). Note that the Dyson equation in 
this form is exact and that many-body approximations 
enter via Sf [G*^]. 



Eq. ( 10 1 is the central equation which will be applied in 



the subsequent Sees. Ill and IV to investigate the effect 
of carrier-carrier correlations in the iV-electron quantum 
dot. However, as the self-energy E*^ appears as a func- 
tional of the Matsubara Green's function G*^(ri,r2; r), a 
self-consistent solution of the Dyson equation is required 
to accurately characterize the equilibrium QD state. The 
corresponding numerical technique is developed in the 
next section. 



III. SIMULATION TECHNIQUE 

In this section, we discuss the computational scheme 
of solving the Dyson equation for the 2D few-electron 
quantum dot specified by Eq. ([3]). Thereby, we proceed 
in two steps: First, we solve Eq. ( [lO| ) at the Hartree- 
Fock (HF) level, see Sec. Ill A[ and, second, we incor- 
porate correlations within the <i>-derivable second or- 



der Born approximation, see Sec. IIIB Throughout, we 
represent G*^ in the r-domain rather than solving the 
Dyson equation in frequency space where G*^(ri, r2;u!) — 
Jgdr G*^(ri, r2; r) e*'^'^ can be obtained by analytic con- 
tinuation, see e.g. Ref.'^Sl. 
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A. Hartree-Fock at zero and finite temperatures 



At mean-field level, the solution of the Dyson equation 



( 10 1 is fully equiv alent to the Hartree-Fock self-consistent 
field methocp2EI] at finite temperatures /?~^. Hence, we 
primarily resort to standard HF techniques and will re- 
cover the uncorrelated Matsubara Green's function, de- 
noted G'°(ri, r), at the end of this section. 

The Hartree-Fock approach leads to an effective one- 
particle description of the QD and gives a first estimate 
of exchange effects. However, as an independent-electron 
approximation, it does not include correlations, i.e. the 
HF total energy is given by E^p = £^cxact - Ecorr- With 
respect to the second quantized Hamiltonian of Eq. Q , 
the effective HF Hamiltonian is obtained by approxi- 
mately replacing the four field operator product entering 
the interaction term by sums over products 4'''' 4' weighted 
by the generalized carrier density matrix p(r,f). This is 
consistent with the mean-field approximation for the two- 
particle Green's function as given in Sec. |nB| and leads 
to 



Hx= //dVd2fV'^(r)[/i°(r)<5(r-r)+S°(r,f)]i^(r) , 



with the Hartree-Fock self-energy 



S",(r,r)= /dV 



, Ap(r',r') 



S{r - f ) 



Ap(r,r) 



(11) 



(12) 



Here, the first (second) term constitutes the Hartree 
(Fock or exchange) contribution. 

Computationally convenient is the introduction of a 
basis representation for the electron field operator ac- 
cording to 

V'(^)(r)=^^W(r)a(^\ ^ e {0,1,2, ...} , (13) 

i 

where the one-particle wave functions (orbitals) fi{v) 

form an orthonormal complete set and a|^^ denotes the 
annihilation (creation) operator of a particle on the level 
i. At this stage, the QD system (111 can be transformed 

with 
being 



into the matrix representation /i^.y = h^j + 
the single particle quantum numbers i and j, hi 
the electron HF total energy, /i^^ the single-particle (ki- 
netic plus confinement) energy and S° -j the electron self- 
energy in mean-field approximation. More precisely, we 
have 



(14) 



d^r(p:(r)(-V^+r^)(p,(r) , (15) 



(16) 



kl 



with the finite (zero) temperature charge carrier density 
matrix Pij{P) = (ajflj) (in the limit /3 — s- oo) in the grand 



canonical ensemble and the two-electron integrals w^j 
defined as 



ij,kl 



d'rd^ 



y^{r — f )2 -I- 



(17) 



Using a ^ 0, the integrals in w^^ can be performed 
analytically in 2D but, in the limit of large anisotropic 
confinement (quasi-lD quantum dot), a truncation pa- 
rameter < Of <C 1 is needed to regularize the (bare) 
Coulomb potential at |r — r| = keeping w^^j j,; finite, 
see e.g. Ref.^°. Alternatively, the parameter a adjusts a 
confining potential in the perpendicular dimension and 
allows (at small rij ) for a transversal spread of the wave 
functiorF^. For the specific choice of the parameter a, 
see Sec. |lVl 

Using standard techniques, we iteratively solve the self- 
consistent Roothaan-Hall equation^22! for the Hartree- 



Fock Hamiltonian hx^ij, Eq. (14 1, 



E 



(18) 



which at finite dimension rif, x n;, (i = 0, 1, . . . , nt — 1) 
yield the numerically exact eigenfunctions (HF orbitals) 
expanded in the form 4>\,i{r) — X]J=o ^ji Vji'"^): Cij S M, 
the corresponding energy spectrum (HF eigenvalues) e° 
and the chemical potential ijP. Consequently, the A^- 
electron quantum dot system is fully characterized by the 
solution (j)\,i(r), e.g. its charge carrier density is given by 



"6-1 



o°(r) = J2 /(/3,e^-/)'/'A,.(r) 



(19) 



j=0 

"6 — 1 



= ^/(Ae°-/)^c,.^,(r), 

i=0 j=0 

where f{l3,e^ — pP) denotes the Fermi-Dirac distribution. 

For numerical implementation of the mean-field prob- 
lem ( [T8| , we have chosen the Cartesian (2D) harmonic 
oscillator states 



</2m,n(r) 



\/2™+" m! n! tt 



Kn{x)nn{y), (20) 



with single-electron quantum numbers i — {m,n), r — 
{x,y) in units of the oscillator length ^q, the Her- 
mite polynomials TLrn{x) and (m -I- l)-fold degener- 



ate energies tr, 



1, m,n e {0,1,2,...}. 



In the ID quantum dot limit, these states reduce to 
the one-dimensional oscillator eigenfunctions tpm{x) = 



(2™to!v^)- 



-1/2, 



As mentioned before, the self-consistent Hartree-Fock 
result generates an uncorrelated Matsubara Green's func- 
tion G'^(ri,r2; r) which yields the same observables. For 
instance, the A^-electron density of Eq. (19 1 is recovered 



from p (r) = G (r,r;0 ) — the energy contributions are 
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discussed separately in Sec. IV Expanding G° in terms 
of the obtained HF basis (/)A,i(r) according to 



Solving the self-consistent Dyson equation 
beyond the Hartree-Fock level 



(21) 



with associated r-dependent real matrix elements g^- (r) , 
leads to the identity 



(22) 



In this subsection, we focus on electron-electron cor- 
relation corrections to the self-consistent Hartree-Fock 
reference stat^^] determined by G°(ri,r2;r). The idea 
is to start from the Dyson equation ( 10 1 in HF orbital 
representation 



which is (band) diagonal only in the HF orbital basis and 
solves the Dyson equation in mean-field approximation 



(23) 



(S° )y = I]° ^■ are defined in correspondence to Eqs. ( 15 ) 
and ( 16 1 with Lpi being replaced by 4>\^i, the Hartree-Fock 
GreeiTs function being denoted as (g*'(T))ij = g^Lij). and 



the charge carrier density matrix due to Eq. ( 16 1 reads 
PijiP) = gij{0~) with notation 0~ denoting the limit 
from below on the contour C. Further, it is (ab)^^ — 
J2k^ikbkj- Note, that in Eq. (17 1 also the two-electron 



integrals are to be transformed into their HF representa- 
tion, and that in the following bold-typed expressions as 
introduced in Eq. ( 23 1 denote matrices with respect to 
the HF basis (t)\,i{r). 



i-dr - hO] g^r) = 6{t) + fdf Sf (r - f) g^^(f) , (24) 

Jo 



Here, the time-independent matrices (h°)ij — ft,? and with the 



(Sf(r)),, 



full, time-dependent Matsubara self-energy 
^ijij) and the equilibrium Green's func- 
tion (g**(T))ij = gfj{T), both obtained by applying the 
orbital expansion of Eq. ( |2ip . An explicit approximate 
expression for SjJ^ including correlation effects is intro- 
duced below, cf. Eqs. (pQpTl. 



First 



we discuss the general solution scheme for 
(24 1. However, we will not consider it in this form. 



Eq. 

Instead, we integrate Eq. (24) inserting Eq. (22 1 and ap- 
plying the anti-periodicity property of g*^(T) . This leads 
to the integral form of the Dyson equation 



\r) 



df df J 
Jo 



f)S^,[g^^](f-f)gnr) = 



]{r) 



g°(r) , 

Sf[g^'^](r)-5(r)Si, 



(25) 
(26) 



where the expression S^(r) according to definition (26) 
implicates the total Matsubara self-energy reduced oy 
the initial (steady-state) mean-field S| — S°[g°(0^)] 
which is not a functional of the full (correlated) Green's 
function gj^^r). In addition, the single-particle energy 
h*^ has already been absorbed in the HF reference state 
g°(r) and thus does not appear explicitly in Eq. (25 1. 



For a more detailed derivation of Eq. ( 25 ) see Appendix 



We highlight, that the integral form of the Dyson 
equation can be parameterized by the second index 
j £ {0, 1, — 1} of the Matsubara Green's func- 
tion gfjir) since the matrix multiplications on the left 
hand side of Eq. ( p5| do not affect this index. Hence, at 
a fixed Matsubara self-energy and discretized r-interval 
[—/3, /3], Eq. (25 1 allows for reinterpretation as a set of Uf, 
independent (but typically large-scale) linear systems of 
the form 



(27) 



^ gfjirp) and (S(-'")),p = 5°-(Tp), respectively. 
The coefficient matrix {A)ip.jq — aij^Tp, Tq) is defined by 
the expression (convolution integral) 

a^J(T,f) = %(5(T-f) (28) 
"1,-1 



in which the integral over r in the Dyson equation (24) 



vanishes due to its replacement by the matrix multiplica- 
tion ylA'(j). In more detail, we need to specify the time- 
discretization of the Matsubara Green's function under- 



where the unknown quantity and the inhomogeneity are 



taken in Eq. (27): First, due to the anti-periodicity prop- 
erty of G*^, we can restrict ourselves to solve Eq. (24) on 
the negative r-interval [— /3,0]. This specific choice orig- 
inates from the fact that in the limit r ^ 0^ the density 
matrix is obtained from ^^{t). Second, the numerical 
treatment must tak e in to account the time-dependence of 
G^^{t). From Eq. ([22| it follows that the Green's func- 
tion is essentially peaked around r = and ±/3. Thus, 
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not an equidistant grid but a uniform power meslP^ 
(UPM) is adequate to represen t the Green's function — 
this method is also used in RefsPSlSSl With a total num- 
ber of n,„ mesh points the dimensionality of the linear 
system AX'-^'> — B'^^^ becomes ntUm x nbn,n- As stated 
above, Eq. (27 1 can only be processed for a fixed self- 
energy \ {t) . This means, that in order to provide 
a self-consistent solution of the Dyson equation, we have 
to iterate the procedure by computing, at each step, a 
new self-energy from the current g^^(T). This loop is 
then repeated until convergence. 

So far, we have no t sp ecified a certain self-energy ap- 



proximation. In Eq. (26 1, one generally can split S*^(t) 



into a mean- field and a correlation part, i.e. 

^[/■^](r) = S{r) S°[/-^(0-)] + ^T"[9'']{r) 



where the Hartree-Fock contribution. 



yO 



kl 



ij,kl 



"il.kj 



(29) 



(30) 



is exact (compare with Eq. ( 16 1) and the correlation part 
Yi'^°"{t), at the second Born level, is given by 



klmnrs 



X3^(r)gr(r)g-(-r) 



M/ 



Here, the first term denotes the direct contribution 
whereas the second one includes the exchange — for de- 
tails see e.g. Refs.^"'^''. Note, that the two-electron in- 
tegrals are given in their Hartree-Fock basis repre- 
sentation. Since the interaction potential M^(l,2) enters 
Eq. (31 1 in second order the present description of charge 



carrier correlations goes beyond the first Born approxi- 
mation of conventional scattering theory. 

When the self-consistency cycle reaches convergence, 
the matrix g^^(T) becomes independent of the initial 
state g°(T) and, in configuration space, the correlated 
Matsubara Green's function of the QD system ^ fol- 
lows as 



G*^(ri,r2;T) 



E 



0L(ri)</>A,,(r2)<?i\'(r), (32) 



where the HF orbitals <j>\j{r) are those obtained in 
Sec. HI A Consequently, correlations are included via 
the r-dependent matrix elements gfjir), t e [— /3, 0], 
which give access to the electron density via p{y) ~ 
G*^(r, r; 0~). We note, that as obtained from the self- 
consistent HF calculation is only one possible reference 
(init ial) state which can be used in the Dyson equation 
(251. Also different types of uncorrelated Green's func- 



tions e.g. obtained from density-functional theory (DFT) 
or orbitals in local density approximation (LDA) are ap- 
plicable if they satisfy the correct boundary conditions. 
For a recent discussion on the relevance for atoms and 
molecules see Ref.'22l. 



iV = 3 


(ID) 














A 


^HP 




-E'2ndB 


Eo 


£hf 




-Eqmc 


13^1 


1 


8.173 


4.621 


8.201 


6.115 


2.321 


-0.235 


7.661 


2 


10.066 


6.124 


10.215 


6.556 


4.405 


-0.747 


9.510 


13^2 


1 


7.043 


4.670 


7.065 


5.027 


2.153 


-0.115 


6.761 


2 


8.790 


6.169 


8.941 


5.311 


3.936 


-0.306 


8.603 


4 


11.732 


8.852 


11.918 


5.920 


6.303 


-0.304 


11.721 


6 


14.387 


11.231 


14.374 


6.712 


7.822 


-0.160 


14.403 


8 


16.790 


13.362 


16.747 


7.514 


9.354 


-0.120 


16.809 


10 


19.005 


15.294 


18.962 


8.257 


10.800 


-0.095 


19.034 


/3 = 10 














-^QMC 


1 


6.615 


4.673 


6.591 


4.645 


1.987 


—0.042 


6.529 


2 


8.480 


6.173 


8.421 


4.966 


3.560 


-0.105 


8.371 


4 


11.667 


8.853 


11.578 


5.817 


5.917 


-0.156 


11.484 


6 


14.374 


11.243 


14.292 


6.710 


7.720 


-0.137 


14.161 


8 


16.787 


13.376 


16.721 


7.534 


9.296 


-0.110 


16.570 


10 


19.004 


15.298 


18.950 


8.285 


10.752 


-0.087 


18.791 


iV = 6 


(ID) 














A 






-E'2ndB 


Eo 


-Ehp 


EcOYY 


-Eqmc 


/3 = 10 


(GS) 














1 


27.600 


9.263 


27.519 


18.613 


9.028 


-0.123 




2 


36.145 


12.195 


35.919 


19.976 


16.289 


-0.346 




4 


50.960 


16.110 


50.440 


23.666 


27.384 


-0.609 





TABLE I: Different energy contributions in dependence on the 
coupling parameter A for the 'ground state' (GS, l3 = 10) and 
equilibrium states {(3 — 2 and 1) of A'^ = 3 and 6 spin-polarized 
electrons in a quasi- ID quantum dot. fi'^ gives the chemical 
potential as obtained from the H artree-Fock calculation with 
total energy -Ejjfi see Sec. 



Ill A 



-B2ndB, Eq,Eb_f, and Ec. 



are computed from the correlated Green's function g (r). All 
energies are in units of Eq and the underlined values pertain 
to the results shown in Figs. |3|5l and Fig. [6] For comparison, 
Eqmc denote the total energy obtained from quantum Monte 
Carlo (QMC) simulations, see also Fig. [5] 



In summary, the presented procedure is valid for arbi- 
trary temperatures l3~^ and arbitrary coupling parame- 
ters A. Thereby, the scope of numerical complexity is de- 
termined by the parameters nt (matrix dimension associ- 
ated with the HF basis size) and rim (time-discretization 
on the UPM) which must be chosen with respect to 
convergence of the QD observables. Corresponding to 
Eq. ( 25 1 it has been found that particularly the particle 
number TV = X]"=o 9iii^~) ^^'^ correlation energy 
i?corr sensitively depend on Um, cf. Eq. (34 1 in Sec. IV 



IV. NUMERICAL RESULTS 

In this section, we report on the numerical results for 
the few-electron quantum dots with = 2, 3 and 6 
charge carriers. At that, we mainly focus on the en- 
ergies and the (accumulated) single-carrier density and 
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I ^ \ I \ I I I L 

1 2 4 6 8 



coupling parameter A 




FIG. 1: (color online) The six energetically lowest HF or- 
bital energies 6° in dependence on A for the three-electron 
QD at different temperatures /3~^. The gray area indicates 
the ground state HOMO-LUMU gap between the occupied 
and unoccupied states. The chemical potential pi" (A) (double- 
dotted-dashed curves) is situated within this gray area. Inset: 
A-depencence of the correlation energy i?corr , see Table |l] 



compare the influence of HF and second Born type self- 
energies, i.e. Eq. (30 1 versus Eqs. (30 1 plus (31 1. 



The energies that contribute to the total energy of the 
QD system are, in addition to the single-particle (kinetic 
[t°] and confinement [v"]) energy Eg = Trh*'g*^(0^) = 
Tr (t"-|-v°) g^(0"_K the mean-field Hartree-Fock and the 
correlation energy^D defined as 



i?HF = xTrSOg*^(0-), 



(33) 




4 6 8 

coupling parameter A 



10 



FIG. 2: (color online) Total energies in dependence on A and 
l3 as given in Table |l] for the three-electron QD. Comparison 
of the [grand canonical] Green's function result (at HF and 
second Born level) with quantum Monte Carlo [canonical]. 
For A = 0, the total energy can be analytically obtained from 
the (grand) canonical partition functio n of t he noninteracting 
system according to standard formulaPES. 



eral, it is 



n,(iV;A,/3)=5-/(0-), 



(35) 



which, in the case of vanishing correlations (G*'^ G^), 
is just the Fermi-Dirac dis trib ution, i.e. = 5^i(0") 



/o, cf. Eq. (221 



drTrSr^(-r)g*^(r) 



(34) 



The total energy is then ii'2ndB — Eq + E^-p + -E-corr- 
For comparison, the total energy with respect to the HF 
Green's function G"(ri,r2;T) will be denoted by -E^p. 
For evaluation of the two-electron integrals needed in 
Eqs. (33 1 and (34 1, we have chosen the truncation pa- 



rameter a = 0.1 in ID and a = in 2D [no divergence of 
the integrals Wij^ki, see Eq. (17 1 in Sec. HI A . 



Moreover, the HF orbital-resolved energy distribution 
functions (level occupation probabilities) ni{N; X, f3) are 
analyzed with respect to correlation induced scattering 
processes of particles into different energy levels. In gen- 



A. Limit of large anisotropy (ID) 

When in Hamiltonian ([T]) the isotropic confinement of 
frequency Uq is replaced by an anisotropic entrapment 
according to Uy^o ^ ujx.Oi the QD charge carriers move 
effectively in one dimension. With the finite regulariza- 
tion parameter a = 0.1 we thereby allow for a small 
transversal extension (perpendicular to the x-axis). That 
is why we will call the system in this regime quasi-one- 
dimensional (quasi- ID). 

In the following, let us first consider the ID version 
of quantum dot lithium'^^ {N = 3) and, hereafter, a QD 
with = 6 confined electrons. For the respective HF 
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FIG. 3: (color online) Thermodynamic properties of = 3 
electrons in the quasi-lD QD at /3 = 10 and A = 2.0. a) Ideal 
and HF energy distribution functions ff, b) change Ui — of 
the HF distribution due to correlations (in percent), and c) 
charge density profiles p(x). The ideal (dotted line) and the 
HF result (solid line) are displayed together with the second 
Born approximation (dashed line). 



and second Born calculations we throughout have used 
ni, = 30 oscillator functions, see Sec. |III A| and the num- 
ber of mesh points njn(u,p) for the r-intcrval [— /3, 0], 
discretized in Sec. |IIIB| was varied between 60 to 140 
in order to achieve convergence and preservation of par- 
ticle number in the Dyson equation (251. Table |l] gives 
an overview of the relevant energies obtained at differ- 
ent coupling parameters and temperatures. Also, we in- 
cluded reference data from quantum Monte Carlo (QMC) 
simulation^2ll for the three-electron QD. 



The low-energetic discrete orbital energies con- 
tributing to the HF reference state g°(T) are shown in 
Fig. [l] in dependence on A and f3. For the quasi ground 
state (GS), /3 = 10, the occupied states z < 3 are ener- 
getically well separated from the unoccupied states z > 3 
by the HOMO-LUMO gap — the energy gap between the 
highest occupied (molecular or Hartree-Fock) orbital and 
the lowest unoccupied (molecular) orbital), see the gray 
area. For temperatures /? < 10, this gap is reduced emi- 
nently for moderate coupling around A w 3 (see the dot- 
ted and dashed lines), while, for A ^ oo, the curves con- 
verge due to the strength of carrier-carrier interactions 
exceeding the incluence of thermal fluctuations. More- 
over, the HF chemical potential /i°(A), situated within 
the HOMO-LUMO gap, is only slightly affected by (3, 
compare the values in Table |T] 

If we now include correlation effects, the HF spec- 
tra {e°} become renormalized — for the discussion see 
Sec. |IV C[ On the QD total energy, the influence of cor- 
relations is as follows: The correlation contribution Emrr 
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FIG. 4: (color online) Same as Fig. [3j but for temperature 
f3 = 2. 
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FIG. 5: (color online) Same as Fig. |3j but for temperature 
13 = 1. 



is negative and increases with temperature but it is non- 
monotonic with regard to the coupling parameter A, see 
Table [l] and inset in Fig. [l] More precisely, the correla- 
tion effects are dominant at moderate coupling, between 
A = 2 and 6, leading to lower total energies i?2ndB (com- 
pared to E'^-p) temperature /? = 10 and to an energy 
increase for (3 — 2 and 1 . 

In general, both, the HF and second Born total ener- 
gies well approach the corresponding [exact] QMC data 
independent of A and (3, see Table |lj In order to properly 
compare the different approximations, the total energy is 
shown in Fig. [2] relatively to the mean-field chemical po- 
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FIG. 6: (color online) Thermodynamic properties of = 6 
charge carriers confined in the quasi- ID QD at = 10 and 
A = 2.0. a) Ideal and HF energy distribution functions f", 
b) change rii — /f of the HF distribution due to correlations 
(in percent), and c) charge density profiles p{x). Labeling as 
in Fig.[3l5l 



tential /i*^. In the whole considered A-regime, the quasi 
ground state energy E^p (dashed curve for P — 10) is 
downshifted by second Born corrections (solid curve) to- 
wards the energies obtained by QMC (dotted curve). At 
A = 2, we identify the best agreement of the correlated 
result with i^QMC- In particular, for A — > 0, the three 
different total energies converge to the value E — f/' = | 
of the ideal QD. At larger temperatures P — 2 and 1, 
the correlated Green's function g*^(T) leads, at moder- 
ate coupling around A w 3, to total energies that are 
significantly larger than the corresponding HF energies. 
This is consistent with the larger absolute values of the 
correlation energy -Econ- given in Table [l] For stronger 
coupling A > 6 (at /3 = 2), £^2ndB then crosses the HF 
value in order to converge to the respective GS curve. 
The comparison with QMC is difficult at finite temper- 
atures: We point out that (already at A = 0) there is a 
general shift due to the usage of different ensemble av- 
erages. Whereas QMC uses a canonical approach, the 
Green's function results emerge from a grand canonical 
picture. Nevertheless, close to the GS, Fig. [2] reveals a 
quite similar behavior in dependence on A. 

From Figs. 3][5 a), one gathers how the QD mean- field 
renormalizes the ideal equidistant energy spectrum 
Ci — i + ^ of the noninteracting system (A = 0), see 
the shifted HF energies e° (open circles) which exactly 
follow a Fermi-Dirac distribution according to Eq. ( 22 1 . 
Correlations due to S™'"' now modify this statistics as 
can be seen from the quantity rii — ff in Figs. [3]|5]b), 
which measures the HF orbital-resolved deviation from 
the Fermi-Dirac distribution (in percent) and shows that 



charge carriers around are being scattered into higher 
HF orbitals (black circles). At larger temperatures, see 
e.g. Fig. |4]b), the change in the occupation probabifity 
rii exceeds 2 % for A = 2.0. Moreover, Pauli blocking in- 
hibits energetically low lying electrons to essentially take 
part in the scattering process — consequently, — is 
small for e° <C /i". 

Figs.pllS^c) visualize the inhomogeneous one-particle 
density p[x) in the three-electron quantum dot (A = 2.0) 
in HF and second Born approximation. For comparison, 
we included also the density of the respective ideal sys- 
tem (blue dotted lines). Being symmetric around a: = 0, 
the density at low temperatures (Fig.[3]and|4]c)) is three- 
fold modulated, and due to the electron-electron interac- 
tions the modulation in p(x) is more intense compared 
to the ideal QD, where the oscillations originate from the 
Pauli principle only. Notably at /3 = 2, correlations sub- 
stantially weaken the density modulation and hence are 
important. At temperature (3 = 1, the correlation effects 
are still present (see the quantity iii — f^) and reveal 
a more smooth, almost monotonic decay of the density 
with a; — > oo. 

For the quasi 'ground state' {(3 = 10) of the QD with 
six electrons, see Table [l] and Fig. [6] respectively, we ob- 
serve similar properties as for the example of quantum 
dot lithium. Here, p{x) becomes six-fold modulated and 
the quantum dot state can be interpreted as a Wigner 
chain of six aligned charge carriers held together by the 
parabolic confinement. However, we note that no or- 
bitals at A = 2.0 are degenerate, and that there is strong 
overlap of the single-carrier (HF) wave functions (px^i^x). 
In contrast to the TV = 3 quantum dot, the influence 
of correlations on the equilibrium state G^^ (see again 
rii — fi) is stronger leading to considerable lowering of 
the total energy. This is consistent with an increased 
number of carriers and hence increased electron-electron 
collision probabilities, compare Figs. |3]and[6]b). 



B. Isotropic quantum dot (2D) 

With no extra restriction on the carrier motion, Hamil- 
tonian ([T|) describes a 2D QD with isotropic parabolic 
conflnement. Here, we report on the obtained Matsub- 
ara Green's function results for the special case of = 2 
electrons, i.e. for spin-polarized quantum dot helium. 
Thereby, we restrict ourselves to very low temperatures 
in order to compare with ground state data available in 
the literature!^. 

For the HF and second Born calculations, the results of 
which are shown in Table [H] the inverse temperature was 
set to (3 = 50, and we used up to = 40 of the energet- 
ically lowest Cartesian oscillator functions, see Eq. (20 1 
in Sec. HI A Further, the uniform power mesh [n 
was chosen as 



IV A 



including essentially more 



in Sec. 
than 100 grid points. 

First, the unrestricted HF energies E^p in Table |ll] 
exactly agree — to more than three decimal places — with 



10 



N = 2 


(2D) 






















X 


-E'cxact 






-B2ndB 


Eo 


-Ehf 


-E'corr 


Kf [%] 


AandB [%] 


C [%] 


pfJ=5 
^QMC 


1 




3.604 


2.885 


3.597 


3.023 


0.584 


-0.010 








3.591 


2 


4.142 


4.168 


3.393 


4.147 


3.082 


1.092 


-0.026 


0.638 


0.121 


80.8 


4.148 


4 


5.119 


5.189 


4.304 


5.117 


3.285 


1.901 


-0.069 


1.367 


0.039 


102.9 


5.123 



TABLE II: Dilferent energy contributions in dependence on coupling parameter A for the ground state (/3 = 50) of jV = 2 
spin-polarized electrons in an isotropic 2D quantum dot. The exact energies for A = 2 and 4 are quoted from Refs.l^^l^si g^^^ 
arise by an exact diagonalization method. -Eqmc gives quantum Monte Carlo reference data computed at the temperature 
13 = 5. Aoo ^ [Ex — -EcxactlZ-Ecxact gives the relative error in %. ^ = (-Erf — -E2ndB)/(-BHF ^ -Ecxact) measures the correlation 
induced improvement of the total e nerg y. All values are given in units of Eq = Hujq. In particular, all three decimal places of 
the HF energies E^p agree with Ref.^. 



the total energies computed in analogous manner by 
B. Reusch et al., Ref.ES. This is a clear indication for 
the HF basis to be large enough. Second, applying the 
second Born approximation we are able to (essentially) 
improve these ground state results. The values of -E2ndB, 
thereby, come quite close to the exact energies, which are 
obtained by numerical diagonalizatioiF^, and the data 
are also in good agreement with the QMC results. 

For instance, for coupling parameter A = 2.0, the in- 
clusion of correlations reduces the relative error A^^ by a 
factor of five — for definition of Ahf and A2ndB see cap- 
tion of Table |llj Hence, with respect to the HF solution, 
this means an improvement of the ground state total en- 
ergy by about ^ = {E^p- E2ndB)/{E^F- Ec^^^a) « 80 %. 
For the strongly correlated case A = 4.0, the calculation 
slightly over-estimates the influence of correlations and 
leads to a total energy lower than the exact value. This, 
most probably, can be improved by increasing the num- 
ber of mesh points for the r-interval which is crucial for 
the convergence of the correlation energy, see integral in 
Eq. (34 1. Nevertheless, the exact ground state energy is 



considerably well approximated. 

As a final remark, we note that our procedure does 
not restore the rotational symmetry of the Hamiltonian 
([ij into the two-dimensional solution. This is due the 
fact that we start from a symmetry broken*^ (HF) ini- 
tial state when solving the Dyson equation to self- 
consistency. For the problem of restoring the symmetry, 
see e.g. Ref.!^. However, the presented results directly 
apply to QDs where impurities naturally avoid the ro- 
tational symmetry and hence lead to symmetry broken 
electron states. 



C. Spectral function 



cate to extract the correlated single-particle energy spec- 
trum as its computation besides Fourier transformation 
usually involves Fade approximation^^. The direct time- 
propagation of the equilibrium state G^(ri,r2;r) solv- 
ing the Keldysh/Kadanoff-Baym equations ^ for the 
two-time NEGF G(l,2) = 0{ti - i2)G>(l,2) - e{t2 - 
ii)G^(l,2) initiates a more systematic approach. Here, 



-1.0 K'^' 
1.0 



1.0 






I ^ 75=0.281 1 


1.0 













1.0 




7;4 = 0.284^ , , , ^ 












1.0 






II 












L ,70=0.181 , 1 







P = 10 



A = 2.0 



-2 2 



6 8 



In HF approximation, the single-particle energy spec- 
trum of the QD consists of discrete levels, see e.g. {e°} 
for the three-carrier system considered in Fig. [T] When 
correlations are included the spectra generally turn into 
continuous functions of energy due to electron-electron 
scattering and provide additional informations such as 
finite line widths or temperature broadening. However, 
from the Matsubara Green's function, Eq. (32 1, it is intri- 



FIG. 7: (color online) Energetically lowest intraband spectral 
functions at (r) (circles) for the one-dimensional quantum dot 
with N = 3 charge carriers at /3 = 10 and A — 2.0. The solid 
lines are the fitting results of the IHC model, Eq. (37 1, with 
exponential damping constants 7^ = rjkVk (solid exponential 
curves). The undamped oscillations (thin dashed curves in the 
background) correspond to the Hartree-Fock approximation, 
where no collisional broadening and thus no damping of (r) 
occurs. 
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FIG. 8: (color online) Spectral function a{u)) (dashed curves filled gray) accumulated from the orbital-resolved functions afe(tj) 
(thin dotted curves) for the quasi one-dimensional QD -with TV = 3 electrons at different temperatures (3~^ as indicated — a) to 
c) A = 1.0, d) and e) A = 2.0. Note that the j/-axis in b), c) and e) is stretched by factor 8. The vertical solid lin es d enote 
the Hartree-Fock energies e°, the dotted lines correspond to the maxima of the inverse hyperbolic cosines ak{uj), Eq. (37 1. The 
numbers at the peak profiles in a) and d) denote the shift of the maxima in ak=2,3{i^) with respect to the HF energies. The 
triangles on the abscissas give the position of the chemical potential f/'. Moreover, the dotted-dashed (double-dotted-dashed) 
curve show the (inverse) spectral weight — for more detailed information see text. 



the time-dependency of G(l,2), which now extends also 
along the real part of the contour C, provides access 
to the electron energies w via the spectral function^^ 
a{uj) = /J*^ dre'"'^a(T), where here r is the difference 
of the two real time arguments in G. 

The orbital-resolved carrier spectral function a(T) is 
given by 

a(r) = z{g> {T-T,T + T)-g<iT~T,T + T)} , (36) 

where T > is a specific point on the diagonal of the two- 
time plane V = [0, oo] x [0, oo] , r = -ts G hV^T, y^T] 
denotes the relative time, and g^{ti,t2) are the contour- 
ordered correlation functions with respect to the HF ba- 
sis. Further, we identify the diagonal (offdiagonal) el- 
ements of matrix a(r) with the intraband (interband) 
spectral functions. 

Well documented computational specifications for 
solving Eq. ^ on V, with in itial condition g^^(T), are 
provided e.g. by Refs-f^f^llMl, We will not give a de- 
tailed description here. Nevertheless, it is instructive 
to note that, as the one-particle energy h" includes no 
time-dependency, the dynamics of G(l,2) can be ob- 
tained by propagation on the ii_2 axes only [instead on 



whole P] or by using the retarded Green's functiorP 
G«(l,2) = i([^(l),^t(2)]+)0(t,_t2). 

As a result it turns out that the QD spectral function 
is not of Breit-Wigner type, i.e. does not obey a distri- 
bution function akk{^) = o.k{^^) oc as it fol- 
lows from a quasi-particle (local approximation) ansatd^ 
with afc(T) = e"^'''^e~''"^, single-particle energy and phe- 
nomenological damping 7. In contrast, the spectral func- 
tion shows clear non-Lorentzian behavior, cf. the circles 
in Fig. [7] which shows the energetically lowest intraband 
spectral functions of the three-electron quantum dot dis- 
cussed in Sec. |IV A[ To this end, at large T 1, we 
have adapted the computed spectral function a^^T) to 
an inverse hyperbolic cosine model^ (IHC) 



afe(r) 



cosh'"' (j/fcT) 



(37) 



which has been demonstrated to yield good results for 
Coulomb quantum kinetics, see Ref.lSSI. The ansatz (l37| 



leaves open a set of three parameters {wfc, 77^, i//c} (_ob- 
tained by fitting), and, in accordance to the numerical 
data, (i) ensures zero slope of 9^eafe(r) at r = and 
(ii) exhibits an exponential decay for large r with a damp- 
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ing constant 7fe = ryfei^fc. The former feature is especially 
missing in the quasi-particle picture. The solid curves 
in Fig. [7] exemplify the good quality of the IHC model 
including properties (i) and (ii) and justify its usage. 

In energy space, the 7fc-induced coUisional broadening 
of the peaks in Uk (to) can be shown to be again described 
by an inverse hyperbolic cosine*^. Explicitly, the accu- 
mulated energy spectrum follows from 

nb—l nb—1 /.4-00 

a{uj) = ^ afc(^) - ^ / dre'^'akir) , (38) 

k=0 k=0 

where in Hartree-Fock approximation one recovers 
a{oj) = X^I-lo^ — e°), compare with Fig.]?] 

For the quasi- ID quantum dot filled with iV = 3 
electrons at coupling parameters A = 1 and A = 2, 
Fig.jsjshows the spectral function a(a;) including all low- 
energetic orbitals at different temperatures . The ver- 
tical solid lines indicate the discrete HF spectra {e°} and 
the gray filled dashed curves show aiui) at the second 
Born level, being composed of the intraband functions 
ak{(jj) which, themselves, are represented by the thin 
dotted curves. Additionally, the positions of the max- 
ima in af^iui) are marked by the vertical dotted lines. 
As a general trend, we observe a shift of almost all 
peaks in a{uj) compared to the HF eigenvalues of the 
charge carriers — in particular, the shifting is dominant 
around the chemical potential . Moreover, the energy 
shifts are accompanied by a state and temperature de- 
pendent coUisional broadening (finite lifetime), where, 
at sufficiently low temperatures, the spectral width is 
small and the intraband functions ak{uj) do not over- 
lap. Close to the ground state, less occupied HF states k 
around ijP, typically, show larger broadening than more 
strongly occupied states, see Fig. |8]a) and d). We note 
that the HOMO-LUMO gap, i.e. the energy gap between 
the highest occupied (molecular or Hartree-Fock) orbital 
and the lowest unoccupied (molecular) orbital appear- 
ing at quasi zero temperature (3 = 10, is reduced by 
electron-electron collisions, and is particularly softened 
at larger temperatures, compare with Fig. |8] b) and e) 
where (3 = 2. This clearly affects the optical absorp- 
tion (emission) spectra of the few-electron quantum dot'^. 
Finally, at even higher temperatures /3 < 1, all spec- 
tra a{u!) gradually become smooth functions with no or 
only few distinct peaks around the chemical potential /i". 
Thereby, also low energies, essentially smaller than e^^g, 
can be adopted by the charge carriers. 

In addition, if in Eq. ( 38 1 the intraband functions (uj) 



are weighted by the respective occupation probabilities 
nfc (or their inverse fifc = 1 — rtfc), one turns the energy 

spectrum a{uj) into J2T=^o^ "-fc afc(cj) (or J2T=o^ afc(^))- 
These quantities allow us to determine with which spec- 
tral weight specific electron energies are (are not) realized 
and thus are (are not) measurable in the correlated QD 
state G*^, see the corresponding dotted-dashed (double- 
dotted-dashed) curves in Fig.jsja) to e). 

Beyond the spectral information, the time-propagation 



of G(l, 2) also allows us to keep track of the accuracy of 
the correlated initial state {t) [solution of the Dyson 
equation (pS])]. Qualitatively, the accuracy can be ex- 
tracted from the temporal evolution of the electron cor- 
relation energji^ which is given by 

EcoAt) = \ jd\l<{v,v-t)-E% , (39) 
/<(ri,r2;ii) = - / d3 W^(l - 3) Gi2(13; 23+)L , 

where G12 is used in second Born approximation. When 
the iterative procedure discussed in Sec. [Ill B| has led to 
convergence and thus to a self-consistent solution g*^(T) 
of the QD Dyson equation, the correlation energy -Ecorr(^) 
must stay constant in time. Consequently, the amplitude 
of any small oscillatory behavior of the correlation energy 
(obtained by propagation) serves as a reasonable estima- 
tor for the error Ai?coir- Throughout, with typical errors 
of less than 5%, such a test has been found to be very 
sensitive and useful to verify the presented results. 

V. CONCLUSION 

In this paper, we have applied the method of 
nonequilibrium Green's functions to study inhomoge- 
neous strongly correlated quantum few-body systems: 
quantum dots with up to six spin-polarized electrons 
in thermodynamic equilibrium. At various interaction 
strengths, the self-consistent solution of the Dyson equa- 
tion at the level of the second Born approximation has 
enabled us to focus particularly on correlation phenom- 
ena. 

Close to the ground state as well as at finite tempera- 
tures, the Born approximation results yield considerable 
improvements for the total energies, the one-electron den- 
sity, and the orbital-resolved distribution functions which 
give access to the electron-electron scattering processes 
being present in the correlated equilibrium state. Fi- 



nally, the discussion of the spectral function in Sec. IV C 
has implied strong influence of correlations on the op- 
tical emission and absorption spectra of the considered 
few-electron QDs. 

Of course, the second Born approximation is a very 
simple model. It neglects both higher order correlations 
(beyond second order in the interaction) and dynamical 
screening (e.g. GW approximation). Nevertheless, com- 
parison of the Born approximation results to first princi- 
ple quantum Monte Carlo and exact diagonalization data 
suggests that this approximation is well capable to accu- 
rately describe the present system. 

Further, the methods discussed in this paper should al- 
low us to study system sizes up to = 12 charge-carriers 
in ID and A'^ = 6 in 2D in a reasonable computer time on 
a single PC. However, the main restriction (limiting fac- 
tor) is basically not the particle number itself but rather 
the number of basis functions, which — together with the 
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discretized r-grid — sets up the large dimensionality of the 
matrices to be computed and processed. At the mean- 
field level, the evaluation of the two-electron integrals, 
Eq. ( 17 1 , and the transformation of which into the HF 
basis are the most time-consuming parts, while solving 
the corresponding Dyson equation is relatively simple. 
For the second Born case, apart from solving the large- 
scale linear system (27) in each iteration, in particular 



in its integral form which can be derived from Eq. (24 1 



corr 
A 



(r), Eq. (|34|), and 
p8|, both needed 



the computation of the self-energy S 
the convolution integrals aij (r, f) , Eq 
with adequate/high accuracy, account for the complexity 
of the calculation. 

Finally and most importantly, the use of NEGFs pro- 
vides a very general approach as one is capable of com- 
puting also time-dependent observables solving the KBE 
(Sec. IIB) for the two-time Green's function G(l,2) un- 
der nonequilibrium situations. Hence, the presented ap- 
proach allows for the extension to other systems such as 
'QD molecules' (assemblies of single QDs) with interdot- 
coupling and time-dependent carrier transport and QDs 
coupled to electronic leads or external (optical) laser field 
sources. 
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APPENDIX: DYSON EQUATION IN INTEGRAL 
FORM 



In order to develop a numerical solution procedure. 



in the following way: Multiplying Eq. (24 1 by the (band 
diagonal) Hartree-Fock state gP{f — rJT^Eq. (22 1, and 
integrating over t G [0, /?] leads to 

[^dT{[-dr - W>]g''{T)}g"{f - t) (A.l) 
JO 

= g°(f ) + fdrfdr {Sf (r - f) g*^(f)} g°(f - r) . 
Jo JO 

Using the identity -dr{g'^' (t) g,°{f - t)} = -drg'^' (t) 
■ gO(f - r) + ^(f _^)g"(f - r) • g*^(T), the left hand side 
of Eq. (A.l) can be written as 



(A.2) 



+ / dr [-a(,_,)g"(f - r) - hO gO(f - r)] g^r) 

JO 

JO 



drS^g^T-r)g-(r), 



(A.3) 



where the first term in Eq. (A.2 1 vanishes due to the anti- 
periodicity property of g*^(r) and in the second term we 



are allowed to insert the Dyson equation ( 23 1 for the 
HF reference state g'^(f — r). Doing so gives us expres- 
sion (A.3). Finally, equating the r.h.s. of Eq. (A.l) with 
(A.3 1 leads over to 



g*^(f)-g"(f) = /dr/dfgO(f-r){Sf(r-f) 
Jo Jo 

-<5(r-f)S°}gnr), (A.4) 



we, in Sec. IIIB have considered the Dyson equation 



which is equivalent to Eqs. ( 25|26| ) in Sec. IIIB with the 
replacements f ^ r, and the notations = and 
S^(r) = Sf(r)-<5(r)S^ 
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